// PairPotentialFunction.cpp: implementation of the PairPotentialFunction class.
//
//////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>
#include "pair_potential_function.h"
#include "DataHolder.h"
#include "nr_cpp/gamma.h"



//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

PairPotentialFunction::PairPotentialFunction(int i_NInput) : _t(0)
{
  NInput = i_NInput;
  nCall = 0;

  _toPower = 4;
  _B = 100;

  Result.Size(NInput);

  // initiate m_atoms array. Note first atom is assumed fixed. 2nd atom fixed along the x-axis
  // and the 3rd atom fixed in the xy-plane

  if ( i_NInput % 3 != 0 )
  {
    cout << "\n\n i_NInput wrong in pair_potential_function\n\n";
    exit(1);
  }

  m_nAtom = i_NInput/3+2;
  vector<double> zeroAtom(3,0.0);
  for (int i = 0; i < m_nAtom; i++)
    m_atoms.push_back(zeroAtom);

  // populate pair potential coefficient sets

  // first to the power 1

  CoefSet set;
  set.n = 1;
  //set.a.push_back(2.4); set.a.push_back(-9.12E-4), set.a.push_back(-1E-3);
  set.a.push_back(2400); set.a.push_back(-9.12E-1), set.a.push_back(-1.0);
  set.b.push_back(0.783); set.b.push_back(0.106), set.b.push_back(1E-5);
  _allCoef.push_back(set);

  for (int n = 2; n <= _toPower; n++)
  {
    CoefSet s;
    s.n = n;
    //int index = 0;
    for (unsigned int i = 0; i < _allCoef[0].a.size(); i++)
    {
      for (unsigned int j = 0; j < _allCoef[n-2].a.size(); j++)
      {
        s.a.push_back(_allCoef[0].a[i]*_allCoef[n-2].a[j]);  
        s.b.push_back(_allCoef[0].b[i]+_allCoef[n-2].b[j]); 
      }
    }
    _allCoef.push_back(s);
  }
}


PairPotentialFunction::~PairPotentialFunction()
{
}


void PairPotentialFunction::updateAtomArray(double *Input)
{
  // update atom array
  m_atoms[1][0] = Input[1];  // assign x value of 2nd atom
  m_atoms[2][0] = Input[2];  // assign x value of 3nd atom
  m_atoms[2][1] = Input[3];  // assign y value of 3nd atom
  int index = 1;
  for (int i = 3; i < m_nAtom; i++)
  {
    m_atoms[i][0] = Input[3*(i-2)+1];
    m_atoms[i][1] = Input[3*(i-2)+2];
    m_atoms[i][2] = Input[3*(i-2)+3];
    index++;
  }
}

/** return value for one two atom interaction term
 *
 * @param n The power of a single pair potential 
 */
double PairPotentialFunction::pairValue(const vector<double>& a1, const vector<double>& a2, int n)
{
  double retVal = 0.0;

  if ( n == 0)
    return 1;

  double fac;
  for (unsigned int i = 0; i < _allCoef[n-1].a.size(); i++)
  {
    fac = 1/(1+8*_allCoef[n-1].b[i]*_t);
    retVal += _allCoef[n-1].a[i]*sqrt(fac)*fac*exp( - _allCoef[n-1].b[i]*fac * (
      (a1[0]-a2[0])*(a1[0]-a2[0]) + (a1[1]-a2[1])*(a1[1]-a2[1]) + (a1[2]-a2[2])*(a1[2]-a2[2]) ));
  }

  return retVal;
}

/** return the part of the derivative of a pair potential which is not individual coordinate dependent
 *
 * @param n The power of a single pair potential 
 */
double PairPotentialFunction::derivOfPartOfPairP(const vector<double>& a1, const vector<double>& a2, int n)
{
  double retVal = 0.0;

  if ( n == 0)
    return 0.0;

        for (unsigned int K = 0; K < _allCoef[n-1].a.size(); K++)
        {
          // the derivative with respect to x, y and z for each atom has
          // the following common part
          double fac = 1/(1+8*_allCoef[n-1].b[K]*_t);
          retVal += 2*_allCoef[n-1].a[K]*sqrt(fac)*fac*_allCoef[n-1].b[K]*fac*exp(-_allCoef[n-1].b[K]*fac*( 
            (a1[0]-a2[0])*(a1[0]-a2[0]) + (a1[1]-a2[1])*(a1[1]-a2[1]) + (a1[2]-a2[2])*(a1[2]-a2[2]) ));
        }
  return retVal;
}



//////////////////////////////////////////////////////////////////////
// functions. Here Input goes from 1 to NInput
//////////////////////////////////////////////////////////////////////

double PairPotentialFunction::Value(double *Input)
{
  // first check that point is not beyong any boundaries
  //checkBoundary(Input);

  double retVal = 0.0;

  // update atom array
  updateAtomArray(Input);

  for (int i = 0; i < m_nAtom; i++)
  {
    for (int j = i+1; j < m_nAtom; j++)
    {
      for (int k = 0; k <= _toPower; k++)
      {
        double fis = bico(_toPower,k);
        double fis2 = pairValue(m_atoms[i], m_atoms[j], _toPower-k);
        retVal += bico(_toPower,k)*pow(-_B, k)*pairValue(m_atoms[i], m_atoms[j], _toPower-k);
      }
    }
  }

  retVal *= pow(-1.0,_toPower+1);

  Result.Input(retVal, Input);

  nCall++;

  return retVal; 
}

// calculates/builds the function value
//
double PairPotentialFunction::Value(vector<double> &Input)
{

  double retVal = 0.0;

  // update atom array
  m_atoms[1][0] = Input[0];  // assign x value of 2nd atom
  m_atoms[2][0] = Input[1];  // assign x value of 3nd atom
  m_atoms[2][1] = Input[2];  // assign y value of 3nd atom
  int index = 1;
  for (int i = 3; i < m_nAtom; i++)
  {
    m_atoms[3][0] = Input[3*i+0];
    m_atoms[3][1] = Input[3*i+1];
    m_atoms[3][2] = Input[3*i+2];
    index++;
  }

  for (int i = 0; i < m_nAtom; i++)
  {
    for (int j = i+1; j < m_nAtom; j++)
    {
      for (int k = 0; k <= _toPower; k++)
      {
        retVal += bico(_toPower,k)*pow(-_B, k)*pairValue(m_atoms[i], m_atoms[j], _toPower-k);
      }
    }
  }

  retVal *= pow(-1.0,_toPower+1);

  Result.Input(retVal, Input);

  nCall++;

  return retVal; 
}

// calculates/builds the derivative
//
void PairPotentialFunction::BuildDerivative(double *Input, double *derivative)
{
  // update atom array
  updateAtomArray(Input);

  for(int i=1; i<=NInput; i++)
  {
    derivative[i] = 0.0;
  }

  // take the derivative with respect to the p'th atom
  for (int p = 1; p < m_nAtom; p++)
  {
    for (int j = 0; j < m_nAtom; j++)
    {
      if ( j != p)
      {
        double factor = 0.0;

        for (int k = 0; k <= _toPower; k++)
        {
          factor += bico(_toPower,k)*pow(-_B, k)*derivOfPartOfPairP(m_atoms[j], m_atoms[p], _toPower-k);
        }
        factor *= pow(-1.0,_toPower+1);

        if ( p == 1 )  // special case of 2nd atom only derivative with respect to x coordinate
        {
          derivative[1] += factor*(m_atoms[j][0]-m_atoms[p][0]); 
        }
        else if ( p == 2 )  // special case of 3rd atom
        {
          derivative[2] += factor*(m_atoms[j][0]-m_atoms[p][0]); 
          derivative[3] += factor*(m_atoms[j][1]-m_atoms[p][1]); 
        }
        else  // starting from atom 4rd - corresponding to p=3
        {
          derivative[3*(p-2)+1] += factor*(m_atoms[j][0]-m_atoms[p][0]); 
          derivative[3*(p-2)+2] += factor*(m_atoms[j][1]-m_atoms[p][1]); 
          derivative[3*(p-2)+3] += factor*(m_atoms[j][2]-m_atoms[p][2]);
        }
      }
    }
  }
}

// used to test if the function is working properly
void PairPotentialFunction::test_function()
{
  ofstream out("output/test_of_pair_potential.txt");
  
  out << "NUMBER OF DIMENSIONS: " << NInput << "  NUMBER OF atoms: " << m_nAtom << endl
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 1;
  double *mInput;
  double *mDeriv;
  newA(mInput, NInput);
  newA(mDeriv, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(10.0*ran1())-5.0;
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }

    out << setw(16) << Value(mInput);

    BuildDerivative(mInput, mDeriv);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << mDeriv[l];
    }

    out << "\n\n";
    
  }
  out.close();

  out.open("output/test_of_pair_potential.cif");
  printCIF(mInput, Value(mInput), out, "fisse");
  out.close();

  deleteA(mInput);
  deleteA(mDeriv);
}

// used to print the solution (this is the old version, the new version is below)
void PairPotentialFunction::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final PairPotential value = " << funcVal << endl;
  out << "Solution coordinates:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void PairPotentialFunction::print_sol(char filename[])
{
  ofstream out(filename);	

  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final PairPotential value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;

  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}

void PairPotentialFunction::printCIF(double *Value, const double ChiSq, ofstream& out, string cifLabel)
{
  double edge = 500;
  out << "data_" << cifLabel << endl;
  out << "_symmetry_space_group_name_H-M '" << "P1" << "'" << endl;
  out << "_cell_length_a    " << edge << endl;
  out << "_cell_length_b    " << edge << endl;
  out << "_cell_length_c    " << edge << endl;
  out << "_cell_angle_alpha   " << 90 << endl;
  out << "_cell_angle_beta    " << 90 << endl;
  out << "_cell_angle_gamma    " << 90 << endl;
  out << "_cell_formula_units_Z     1" << endl;
  out << "loop_" << endl;
  out << "    _atom_site_label" << endl;
  out << "    _atom_site_fract_x" << endl;
  out << "    _atom_site_fract_y" << endl;
  out << "    _atom_site_fract_z" << endl;
  out << "    _atom_site_type_symbol" << endl;


	  for (int iAtom = 0; iAtom < m_nAtom; iAtom++)
	  {
		  out << setw(4) << "C" << iAtom << " ";
		  out.precision(5);
		  out << setw(10) << m_atoms[iAtom][0]/edge << " ";
		  out << setw(10) << m_atoms[iAtom][1]/edge << " ";
		  out << setw(10) << m_atoms[iAtom][2]/edge << setw(4) << "C" << endl;
	  }


  out << "#END" << endl;
}

// distance from know global minimum which is located at the origin
double PairPotentialFunction::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += solution[i]*solution[i];
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}